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SUMMARY 


A unified formulation is presented based on the full potential framework coupled with an ap- 
propriate structural model to compute steady /unsteady flows over rigid/flexible configurations for 
across the Mach number range (subsonic to supersonic). The unsteady form of the full potential 
equation in conservation form is solved using an implicit scheme maintaining time accuracy through 
internal Newton iterations. A flux biasing procedure based on the unsteady sonic reference con- 
ditions is implemented to compute hyperbolic regions with moving sonic and shock surfaces. The 
wake behind a trailing edge is modeled using a mathematical cut across which the pressure is satis- 
fied to be continuous by solving an appropriate vorticity convection equation. An aeroelastic model 
based on the generalized modal deflection approach interacts with the nonlinear aerodynamics and 
includes both static as well as dynamic structural analyses capability. Figure 1 shows a schematic 
of the coupling process. Some of the topics to be discussed are 1) mechanism for coupling the 
aerodynamic and the aeroelastic equations, 2) update of geometry and grid at each time level to 
maintain time accuracy, and 3) prediction of flutter point dynamic pressure. 

Results are presented for rigid and flexible configurations at different Mach numbers ranging 
from subsonic to supersonic conditions. The dynamic response of a flexible wing below and above 
its flutter point is demonstrated. 



Fig. 1. Schematic of CFD /Aeroelastic Coupling 
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AERODYNAMIC FORMULATION 


The conservation law form of the time-dependent full potential equation cast in an arbitrary 
coordinate system (r, (,»;,£) is written as 


Pr + + Fjf + = 0 


where 




E = (p( t + E( x + F Cv + GG)/J 
F = (pnt + Er) x + Frjy + Grjz)/J 
G = (p( t + E( t +F(, + Gt z )/J 

and J, Jacobian for the transformation,*^ given by 


( 1 ) 


J = d(T,(,ihO/d(t i x t y,z). 

The density p is given by Bernoulli’s law 

P y - , = l- :L p-Ml[24, r + (V + (,) 

+ (V + •/,)*, + (W + (')t ( -l] 

where U, V, and W are the contravariant velocities. 

Associating the subscripts i,j ,k with the directions, a numerical approximation to Eq. (1) 

can be written in the semidiscrete conservation law form given by 

( j) r + (-^«,j,*+i /2 ~ F i,j,k- 1 / 2 ) 

+ (^«,i+i/2,fc - (3) 

+ ( G i+ i/2,j,k ~ G.-i/2,j,*) = 0 

where E,F,G are representative fluxes approximating the real fluxes E,F,G of Eq. (1). 


177 



NEWTON ITERATION 

In terms of the velocity potential, Eq. (3) can be represented as 

W) = 0 (4) 

where <f> is the unknown to be solved at every grid point (i,j,k) in the current (n + 1) time plane. 
The Newton iteration for solution to Eq. (4) is 

* ( « + (fL: A *=° < 5 > 

where <f>* is the currently available <j> at the (n + 1) level and = <f> - <f>*. At convergence, 
will approach zero driving <f> to <j> n+1 . (See Fig. 2). 
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Fig. 2. Update of <t> Based on Newton Iteration 
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TIME LINEARIZATION 


Treatment of ^ ( j) in Eq. (1) 


a <«•-*•){ (5)” + ‘ - (5)” }-»i{(5)“- (5)"“} 

dr \J / a\ Arj — 6i(Ari -f Ar 2 ) 

where 

ai = (An + At 2 ) 2 ; Ari = r n+1 — r n ; Ar 2 = r n - r" -1 . 
The unknown quantity in Eq. (6) is p n+1 . Following Eq. (5), this is written as 


where A<f> — <f> — <f> m , and 


P = P(<t>*) + 



A<f> 




is a differential operator. The speed of sound is denoted by a. 



FLUX BIASING 


Equation (3) requires evaluation of E, F, and G at various spatial half node points. As 
mentioned earlier, E represents E appearing in Eq. (1). The fluxes are defined in the following 


manner. 


~ jj ~ y * „w 

E = P~J ! F = P~J ; G = P~j ( 9 ) 

where p is an upwind biased value of density, designed to produce the necessary artificial viscosity 
for treatment of hyperbolic regions. 


P = ~ p?) ± 


(U d V d " W d \ , 

\c?a( + <?av + q d(j ^ . 


where Q = VU 2 + V 2 + W 2 . 

The quantity ( pq)~ appearing in Eq. (10) is defined to be 

(pq)~ = pq — p*q * if q > q* 

= 0 if q < 9* • 

The quantities p 6 q*< p*, and q* represent sonic values of the flux, density, and total velocity, 
respectively. These quantities are evaluated at half node grid points. (See Fig. 3). 



Fig. 3. Notation for Flux Biasing 
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AEROELASTIC MODEL 


The aeroelastic model is based on the generalized modal approach. In the physical space, the 
structural equation is written as 

[m]z+[c]z+ [*] = {/} 

where m, c, k, and / represent the mass, damping, stiffness, and force, respectively. The structural 
deflection is given by z. Using the generalized mode shape, <f>, one can rewrite the above equation 
as 

"{«} + [<?]{«} + [*]{«}-{*■} 

where M, C, and K are generalized mass, damping, and stiffness matrices and { F } is the gener- 
alized aerodynamic force. These are defined by 

M = <f> T rruf> ; K = <f> T k<f> , 

C — <t> T c<f> , F — <f> T f and z = <f>q , 

where q is the generalized deflection. 
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STATIC FLEXIBLE CASE 


The flexibility of the structure comes into the calculation through stiffness ( K ) and generalized 
mode shapes The steps involved in computing the static flexible equilibrium aerodynamics are 
given below: 

1. Solve rigid aerodynamics for lift force. 

2. Solve Kq = F =» {z} = [<j>] {?}. 

3. Increment {z} in smaller steps and compute modified aerodynamic force, Figure 4. 

4. Repeat Step 2. 

5. Continue until {z} and lift converge. 


z 



Fig. 4. Incremental Shape Change for Static Flexible Calculation 
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DYNAMIC FLEXIBLE CASE 


The aerodynamics and the structural response are computed in a time accurate fashion using 
internal Newton iterations. The various steps involved in this computation are 

1. Compute rigid aerodynamics. 

2. Compute static flexible for a given dynamic pressure Q. 

3. Set up an initial perturbance by perturbing either q or q of any mode. 

4. Solve Mq + Cq + Kq = F =► {z} . 

5. Compute x T ,y r ,z T on the surface. 

6. Define x T ,y T ,z T for each field grid point allowing their values to go to zero at outer boundary, 
Figure 5. 

7. Compute new grid location x n+1 = x n + x r A r, • • • and new metrics. 

8. Converge aerodynamics using internal Newton iteration. 

9. Repeat Step 4. (For prediction of flutter Q gradually vary Q and monitor generalized modal 
deflections.) 



Fig. 5. Grid Update for Dynamic Flexible Calculations 



OTHER FEATURES OF AERODYNAMIC/ AEROELASTIC PROCEDURE 

• Aerodynamics solved by triple approximate factorization scheme ^ — L^L^L V 

• Wake cut modeled by unsteady vorticity convection equation 

• Implicit boundary conditions 

• Internal Newton iteration for time accuracy. 

Code Capabilities 

• Static Rigid 

• Dynamic Rigid 

• Static Flexible 

• Dynamic Flexible 
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FLEXIBLE WING MODEL 


Figure 6 shows the planform shape of a flexible low aspect ratio fighter wing. A sectional 
airfoil profile is also shown in Fig. 6. Note that there are two bumps in the lower section which axe 
due - to the leading edge and trailing edge control surface actuators. This model was designed and 
built by Rockwell and currently is being tested in NASA Langley Research Center’s 16 foot TDT 
wind tunnel. 
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Fig. 6. Planform Shape of a Flexible Fighter Wing 
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NATURAL MODE SHAPES 


Figure 7 shows the first ten natural frequency mode shapes for aeroelastic analysis. Mode 
shapes are prescribed at every grid point on the wing surface. 


ORIGINAL PAGE IS 
OF POOR Q'.'A’ itv 



MODE 1 
6.33 Hz 


MODE 2 MODE 3 

10.27 Hz 13.92 Hz 


MODE 4 MODE 5 

22.73 Hz 32.73 Hz 



MODE 6 MODE 7 MODE 8 MODE 9 MODE 10 

35.48 Hz 40.71 Hz 52.75 Hz 54.03 Hz 55.96 Hz 


Fig. 7. Selected Mode Shapes for Aeroelastic Analyses 
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RESULTS FOR STATIC FLEXIBLE CASE 


Figure 8 shows results for a static flexible computation at Moo = 1.15, a = 6° for two different 
dynamic pressure conditions. Q = 0 corresponds to the rigid case. The deflected shape of the wing 
for Q = 288 psf shows the nose down rotation of the airfoil near the wing tip resulting in tip load 
reduction. This is essential for the wing to be aeroelastically stable. 








Fig. 8. Static Flexible Computation, Moo = 1.15, a = 6° 
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C L VERSUS a FOR STATIC FLEXIBLE 

The Cl versus a curves at two different Mach numbers are shown in Fig. 9. Table 1 shows the 
change in total lift coefficient between rigid and flexible cases for different Mach numbers, angles 
of attack, and Q conditions. 


M 

"1 

a 

<3(psf) 

^te#t 

^kflex 

rigid 

1.05 

1.6 

335 

0.0785 

0.0800 

0.1450 

1.15 

0.5 

230 

0.0433 

0.0400 

0.0685 

0.90 

0.0 

220 

0.0446 

0.0441 

0.0637 


Table 1. Comparison of Lift Coefficients between Rigid and Flexible Calculations 



Fig. 9. Cl versus a for the Flexible Wing 
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TOTAL LIFT COEFFICIENT, C L 


STATIC FLEXIBLE WITH DEFLECTED CONTROLS 


Figure 10 shows a similar static aeroelastic calculation with leading and trailing edge control 
surfaces deflected. Also shown are the surface chordwise pressure distributions at various span 
stations and a cross section of the wing with deflected control surfaces. The span stations where 
the deflected control surfaces are present are drastically altered. 



Fig. 10. Static Aeroelastic Computations with Deflected Control Surfaces 
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DYNAMIC FLEXIBLE CASE 


Figures 11-13 show results for dynamic aeroelastic computations at three different dynamic 
pressures (Q). Figure 11 is for Q < ^Flutter- The structure is aerodynamically stable as shown by 
the decaying amplitudes of various quantities. Figure 12 is exactly at the flutter dynamic pressure. 
The value for the flutter Q was obtained by the nonlinear aerodynamics/ aeroelastic code through 
numerical search. Exactly at the flutter Q, the dynamic response of the structure does not decay as 
shown by results of Fig. 12. Figure 13 illustrates the structural response above the flutter Q. The 
structure is aerodynamically unstable as shown by the growing amplitudes of various quantities. 



Fig. 11. Dynamic Flexible Computation - below Flutter, M 0 0 = 1.15, a = 0°, Q = 360 psf 
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DYNAMIC FLEXIBLE CASE (CONTINUED) 


- 0 . 012 , 

- 0 . 0181 - 




a 



GD ~ GENERALIZED DEFLECTION 



Fig. 12. Dynamic Flexible Computation - at Flutter Point, Mo© = 1.15, a = 0°, Q — 490 psf 



GD ~ GENERALIZED DEFLECTION 
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Fig. 13. Dynamic Flexible Computation - above Flutter, M = 1.15, a = 0°, Q = 550 psf 


191 



